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NOMENCLATURE 


o 

C£ Friction coefficient, r /( 0.5 pM ) 

F Nondimens ional velocity in radial direction (equation 2.11) 

Fr Froude number, V//"g7 

-* “* 2 

g Gravitational acceleration, - 9.81 j m/s 

G Nondimensional velocity in angular direction (equation 2.12) 

o 

h Convective heat transfer coefficient, V/m K, 

= { q v/( T v' T b)’ for heating 
1 q w /( T w~ T sat) ’ for eva P° ration 

-* 

i Unit vector in angular direction 


j 

k 

K 

k 

n 


Nu 

P 

Po 

Pr 

q 


^loss 

r 

R 


Unit vector in the direction normal to the plate 
Cell index in radial direction 
Thermal conductivity, W/mK 

Unit vector in z- direction 

Coordinate normal to the free surface, or number of grid cells in 
radial direction 

Unit vector normal to the free surface 

Nusselt number, = htf/K 

Static pressure, Pa 

Ambient pressure, Pa 

Prandtl number 

Heat flux, V/m^ 

Q 

Volumetric flow rate, m /s 

3 

Local rate of volume loss at a cell next to free surface, m /s 
Radial coordinate, m 
Radius of curvature 
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* t 


Re Reynolds nubmer, M6/u 

t Time, s 

t Unit vector tangential to the free surface 
T Temperature , °C 

u Velocity in angular direction, m/s 

v Velocity in the direction normal to the plate, m/s 

-4 

V Velocity vector, m/s 

w Velocity in the z- direction, m/s 

Vin Average radial velocity at entrance, m/s 

x Coordinate in angular direction, m 

y Coordinate normal to the plate, m 

z Coordinate in the radial flow direction (Fig. 1), m 


Greek Symbols 

2 

a Thermal diffusivity, m /s, or relaxation factor 

6 Film thickness, m 

* 

6 Film height corresponding to Fr = 1, m 

( Nondimens ional coordinate perpendicular to the plate (eqn. 10) 

2 

v Kinematic viscosity, m /s 

3 

p Density, kg/m 

<r Surface tension, N/m 

2 

r Shear stress at wall, N/m 

2 

t Stress tensor, N/m 

hi Angular velocity, rad/s 
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Subscripts 


b Mixed mean (bulk) condition 
c Cartesian components of velocity 
e Free surface 

in Condition at entrance 

out Condition at exit 

r Component in radial direction 

sat saturation condition 
w Condition at solid wall 

<j) Component in angular direction 
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SECTION I 


FLOV OVER A STATIONARY DISK 


1.1 SUMMARY 

Improvements in the theoretical model and computational procedure for 
the prediction of film height and heat- transfer coefficient of the free 
surface flow of a radially- spreading thin liquid film adjacent to a flat 
horizontal surface of finite extent are presented. Flows in the presence 
and absence of gravity are considered. Theoretical results are compared to 
available experimental data with good agreement. In the presence of 
gravity, a hydraulic jump is present, isolating the flow into two regimes: 
supercritical upstream from the jump and subcritical downstream of it. In 
this situation, the effects of surface tension are important near the outer 
edge of the disk where the fluid experiences a free fall. A region of flow 
separation is present just downstream of the jump. In the absence of 
gravity, no hydraulic jump or separated flow region is present. The 
variation of the heat- transfer coefficient for flows in the presence and 
absence of gravity are also presented. 
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1.2 INTRODUCTION 


The radially spreading flow of a thin liquid film is encountered both 
in the case of axisymmetric discharge of liquid through a thin slot and in 
the case of impingement of a liquid jet on a solid surface. These 
configurations are encountered in many engineering devices like spray type 
heat exchangers, cooling towers, spin coating deposition on metal surfaces, 
and rotating condensers and evaporators. The appropriate design of these 
systems require an understanding of the structure of the free surface, the 
liquid flow field and associated friction and heat- transfer coefficients. 
These quantities are usually influenced significantly by the magnitude and 
orientation of the gravitational body force. Therefore, information on 
these flows in the presence and absence of gravity are required for 
successful design and operation of these devices on earth as well as in 
outer space. 

One of the earliest works on the radially spreading flow of a thin 
liquid film was performed by Watson [1] . The film was formed by 
impingement of a liquid jet on a flat horizontal surface. Both analytical 
and experimental studies of laminar and turbulent flows were made. A 
hydraulic jump was found to be present under some flow conditions. The 
analysis covered the supercritical region before the jump where four 
different flow regimes were identified and analyzed using the Pohlhausen 
integral method and a similarity transformation [2] . In the subcritical 
region, the film height was assumed to be constant. For a given location 
of the jump, the subcritical height downstream of the jump could be also 
predicted. 
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The heat transfer to a thin film formed by an impinging liquid jet was 
considered by Chaudhury . [3] , whose analysis basically followed the ideas 
introduced by Watson [1] except that the energy equation was also solved. 
In the region away from the impingement location where a similarity 
solution of the momentum equation was possible, the energy equation was 
solved in a closed form including the effects of viscous dissipation. An 
approximate solution for the simultaneously developing hydrodynamic and 
thermal boundary layers was also obtained using the integral method where 
the temperature profile was approximated by a fourth- order polynomial. 

In a series of studies, Azuma and Hoshino [4] presented the behavior 
of the radially spreading thin film discharging from a cylindrical gap 
formed between a circular nozzle and a flat glass plate. Their 
experimental research determined criteria for laminar- turbulent transition, 
the liquid film thickness, the velocity profile across the film, stability 
criteria, and wall pressure fluctuation. The transition from laminar to 
turbulent flow was established by observing the nature of surface waves 
with a stroboscope. Several types of waves were identified, and it was 
proposed that transition occurs when concentric waves change to 
lattice- shaped waves. This transition depended on the rate of flow, nozzle 
diameter, discharge gap and radial location. The liquid film thickness was 
measured with a needle gage. For all the flow rates considered, the film 
thickness decreased after exit from the nozzle, attained a minimum and 
increased further downstream. The velocity profile across the thickness of 
the film was measured using a laser doppler velocimeter. The pressure 
fluctuation was measured by a pressure transducer which gathered signals 
from two hypodermic needles embedded on the surface of the disk. Most 
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experimental measurements were supported by simple theoretical analyses. 
They identified three regions of the flow: a region just downstream of the 
nozzle where laminar boundary layer develops but has not reached the free 
surface, a region where boundary layer reaches the free surface but flow is 
laminar, and finally a region where flow is turbulent. 

Experimental measurements of local heat transfer coefficient when a 
liquid jet impinges a solid wall and spreads out as thin film was reported 
by Stevens and Vebb [5] . The parameters considered were flow rate, jet 
diameter and nozzle- to- plate spacing. Correlations for local and average 
Nusselt numbers were determined in terms of jet Reynolds number. In a 
later study, Stevens and Vebb [6] investigated the effects of jet 
inclination on the distribution of the local heat- transfer coefficient. 

In the abovementioned studies of Azurna and Hoshino [4] and Stevens and 
Vebb [5,6] , attention was focused only on the supercritical flow, and no 
jump was found to be present. A jump is possible when the flow rate is 
relatively small and disk diameter is relatively large so that transition 
occurs before the film gets washed out of the disk. Vatson [1] and 
Chaudhury [3] identified the presence of a jump, but their analysis again 
focused the supercritical region upstream of the jump. A systematic 
investigation of supercritical and subcritical flows including a hydraulic 
jump, was carried out by Thomas et al. [7,8] and Rahman et al. [9-11]. In 
these studies, the film was generated by discharge of liquid through an 
axisymmetric slot from a pressurized container at the center of the disk. 
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Thomas et al. [7] numerically predicted the film height variation 
using a simple one- dimensional model- Neglecting any variation of velocity 
across the thickness of the film, the continuity and momentum equations 
were integrated across the thin film to develop a single equation for the 
film velocity. Results were obtained for both stationary and rotating 
disks for a number of inlet Reynolds, Froude, and Rossby numbers. In the 
presence of the jump, the outlet Froude number was always assumed to be 
unity to simulate a situation where the fluid experiences a free fall over 
the edge of the disk due to gravity. It was found that the jump moves 
downstream and may get washed away entirely with an increase in the flow 
rate, angular velocity and inlet Froude number. Thomas et al. [8] 
presented experimental data for the liquid film height and a photographic 
study of surface waves for both stationary and rotating disks. The 
measurement confirmed the existence of a hydraulic jump for stationary disk 
which was washed away when the disk was rotated. With the increase of 
rotational speed, the surface waves changed from wavy- laminar pattern to 
radial- wave pattern which carried fluid in rivulets on top of a thin liquid 
sheet adjacent to the disk. 

The studies by Rahman et al. [9-11] solved the flow field numerically 
by a finite- difference method using a boundary- fitted coordinate system. 
No assumptions regarding the velocity variation or friction coefficient 
were required. Since the height of the free surface was dependent on the 
flow conditions and was not known ahead of time, an iterative procedure had 
to be used to determine the correct location of the free surface. The 
original method presented by Rahman et al. [9,11] is termed as the 
’pressure optimization method’ . In this method, the shape of the free 
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surface was assumed to be described by an algebraic equation with two or 
more arbitrary constants. The constants were optimized using an exhaustive 
search technique that minimized the difference between the computed free 
surface and ambient pressures. When a jump existed in the flow field, the 
supercritical and subcritical flows were computed separately, and the 
solutions were matched at the location of the jump. The second method 
presented by Rahman et al. [10] is termed as the ’porous wall method’. In 
this method, the free surface is assumed to be a permeable wall through 
which fluid particles may leave or enter the control volume depending on 
the difference between the fluid and ambient pressures. The shape of the 
surface is corrected in successive iterations until the free surface 
conforms to a streamline and the penetration through the surface reduced to 
zero. This essentially satisfies the kinematic condition on the free 
surface. This method computed the whole flow field as a single domain. 

The purpose of the present investigation is to critically investigate 
the merits and drawbacks of the existing computational methods as well as 
significant improvements for the prediction of thin film flow as discharged 
from a pressurized vessel at the center of a horizontal stationary disk. 
The significant improvements sought are assumptions about boundary 
conditions and inclusion of surface tension at the outlet, an improved 
iterative procedure for determining the film height, and better criteria 
for the estimation of error in computed results. Furthermore, the 
theoretical predictions are also compared with experimental data to 
establish their validity. The details of the flow structure and 
heat- transfer coefficient are also determined for flows in the presence and 
absence of gravity. 
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1.3 EXPERIMENTAL DATA 


The details of the experimental setup and the method of data 
acquisition and its uncertainty estimates are described by Thomas et al. 
[8] . Measurements of film height were made for a thin liquid film 
emanating from a pressurized vessel at the center of a circular horizontal 
disk. The diameter of the disk was 406.4 mm and the diameter of the collar 
through which the discharge was made was 101.6 mm. The gap height between 
the bottom face of the collar and the top of the disk was 0.267 mm, which 
is also the initial thickness of the discharging liquid film. Water at an 
average temperature of 22°C was used as the test fluid. The closed loop 
experimental system maintained an isothermal fluid that discharged 
axisymmetrically on the disk. A noncontact capacitance sensor was used to 
determine the thickness of the liquid film. The probe was moved along a 
radius to determine the distribution of film height. The sensing area had 
a diameter of 11.28 mm, providing reasonably accurate measurements of the 
film height except for the location of the hydraulic jump. The film height 
distribution were determined for 5 different flow rates (7, 9, 11, 13, and 
15 LPM) . Distinct supercritical and subcritical regions isolated by a 
hydraulic jump were present only for the three lower flow rates. For both 
13 LPM and 15 LPM, the subcritical region was washed away from the disk and 
a partial jump was seen near the outlet. In order to avoid any misleading 
conclusion, only data for 7, 9, and 11 LPM will be used for comparison with 
theoretical predictions. Comparison of these experimental data with any 
theoretical prediction was not attempted in any earlier work, and it is 
therefore one of the primary objectives of this paper. 
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1.4 NUMERICAL PREDICTIONS 


The radially spreading flow over a stationary horizontal disk of 
finite extent in the presence and absence of gravity is schematically shown 
in Fig. 1.1. In the presence of gravity, the flow may encounter a 
hydraulic jump where incoming supercritical flow transforms to a 
subcritical flow downstream of the jump. Beyond the disk, the fluid 
encounters a free fall due to gravity. In the absence of gravity, the flow 
remains entirely supercritical and no jump is possible. At the exit from 
the disk, the flow discharges radially. The equations governing the 
conservation of mass, momentum and energy in a thin liquid layer involving 
a Newtonian, constant- property liquid can be written as 


V • V = 0 


DV = . 1 V p + , V 2 V + g 


DT „ u 2 t 
Ut = 0 V T 


The boundary conditions are given by: 


V = V in- T = V *=*i. 


( 1 . 1 ) 


( 1 . 2 ) 

(1.3) 


at entrance: 


(1.4) 



Free surface 



(a) Flow under zero gravity 


Free surface 



(b) Flow under normal gravity 

Fig. 11 Schematic of the flow systems 
and coordinates for computation 


at exit: 


f In the absence of gravity: 
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cE 


= 0, at r = r 


out 


(1.5a) 


In the presence of gravity 
6 = (equation 1.11), at r 

-* = \uv m = °» at y = °» 

-4 

at solid wall: V = 0, T = 


r out 

r out 


< r < r 


e 


(1.5b) 


( 1 . 6 ) 


-4 -4 

at free surface: r • n = 0, r • t = 0, 


d 6 31 _ n 

Hr w * 7E 


(1.7) 


Here n is the coordinate directed normal to the surface under 
consideration. When surface tension is negligible and pressure is the most 
significant normal stress, the r • n condition can be written as p = p Q , 
where ’p 0 ’ is the static pressure of the surrounding medium. 

The boundary condition at the exit depends on the existence of the 
gravitational body force. In the absence of gravity, the flow is entirely 
supercritical and the film height at the outlet does not need to be 
specified. In the presence of gravity, when a hydraulic jump is present in 
the computation domain, a subcritical flow is present in the vicinity of 
the outlet. Then, the outlet film height need to be specified. In the 
study of Thomas et al. [7] as well as Rahman et al. [9,10], it was assumed 
that the Froude number at the outlet (r = r^) was unity, since the 
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subcritical flow becomes critical and finally supercritical as it falls off 

the disk. Therefore, in the absence of any other better approximation, the 

flow at r = r Qut was assumed to be critical. However, in reality, the 

fluid does not experience the fall until it is swept away from the disk as 

shown in Fig. 1.1. Moreover, the experimental measurements of Thomas et 

al. [8] suggests that the flow in the vicinity of the exit is entirely 

subcritical. Table 1.1 shows a comparison between the experimentally 

* 

measured film height and the critical height, 6 predicted by the Fr = 1 

* 

condition at the exit. Note that 6 can be calculated from the continuity 
equation: 


* 

5 


r q f /3 

21 r out 


( 1 . 8 ) 


The measured height corresponded to the data point closest to the edge of 
the disk where measurement was reliable. It appears that the actual film 
height is about 4 times the critical height. Therefore, even though a free 
falling condition exists just downstream of exit, the film is held up by 
some sort of external force. 


It can be recalled here that in all previous computations [7,9,10] , 
the surface tension was neglected in all regions of the flow, which may not 
be correct near the exit since an appreciable curvature of the free surface 
is encountered due to the rapid turning of the flow. In the presence of 
surface tension, the fluid and ambient pressures can be related by the 
equation 
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Table 1.1: Comparison of the Calculated and Measured Film 

Height at the Outlet in the Presence of Gravity 


Critical Height 


Flow rate 

Measured 

Height 

(6 ) from 
eqn. (1.8) 

(LPM) 

(mm) 

(mm) 

7 

4.420 

0.981 

9 

4.318 

1.160 

11 

4.369 

1.337 


Calculated 

Deviation Height (^ ou ^.) Deviation 
of S from eqn. (1.10) ^ou t 


00 

(mm) 

(7.) 

-77.81 

4.377 

-0.95 

-73.15 

4.479 

3.73 

-69.39 

4.582 

4.88 
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( 1 - 9 ) 


Pe- 


P 


<r 

a " E 


where R is the radius of curvature. If the pressure is assumed to be 
hydrostatic in nature, this equation can be expressed in terms of film 
height as 




( 1 . 10 ) 


where ^ is the actual film height at the exit end of the disk (r = 
* 

r . ) , and 6 is the critical height predicted by the Fr = 1 condition at 
that location. To solve equation (1.10), an estimate for R is required, 
since that is unknown in general. The visual observations showed that when 
the flow changes its direction from horizontal to vertical, the radius of 
curvature is of the same order as the film height. Substituting R = 0.5 
and solving for the film height from equation (1.10), we obtain 



* 

ft + 

P - 

* 




(l.n) 


Equation (1 .11) , even though quite approximate in nature, gives a good 
estimate of film height at r = r This is also demonstrated in Table 
1.1, where the height calculated by equation (1.11) is compared with the 
measured height. The deviation is within 5% which is reasonable 
considering the simplicity of formulation. Therefore, whenever a jump was 
present, the film height at the exit end of the disk was computed from 
equation (1.11). The liquid film experiences a free fall as it exits from 
the disk at r = r Based on the experimental observations and 
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comparisons seen in Table 1.1, it was assumed that the free surface forms a 
quarter circular arc with radius R and then extends vertically downward. 
The radial location of the free surface is given by 


r e = r out + °' 5 'out 


The uniform velocity of the free falling liquid was calculated from the 
conservation of mass. 


The coordinate system used in the present formulation is also shown in 
Fig. 1.1. the r-axis is directed along the plate radius and y-axis is 
directed perpendicular to the plate. This r-y system was used for the 
formulation of the one- dimensional computational scheme where the 
continuity and momentum equations were integrated across the y- direction. 
For a complete two-dimensional solution, a body- fitted coordinate system 
[9] was used. In this system, the coordinates were defined along lines 
connecting the centers of the adjacent grid cells. The z-axis lies along 
the main direction of the flow, and the y-axis across the thickness of the 
liquid layer. Comparing the r-y and z-y system, it can be noticed that 
the y-axis remains unaltered, whereas the z- coordinate is related to r by a 
geometrical factor. 

The present one- dimensional computational scheme is similar to that 
developed by Thomas et al. [7] except that the film height at the outlet 
was calculated using equation (1.11) instead of the critical height. In 
this scheme the governing equations for mass and momentum (1-2) were 
integrated assuming a uniform velocity across the thickness of the liquid 
layer, and were combined to give a single equation for film velocity. 
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Since the velocity was explicitly specified at r = r Qut ,the free falling 
flow beyond that location could be ignored. Even though the steady- state 
solution was desired, the transient term was retained in a simplified form 
to devise an explicit marching scheme. The friction coefficient, c^, was 
estimated using a two-part formulation: Blasius- solution for regions where 
the thickness of the boundary layer was smaller than the film thickness, 
and the coefficient obtained from a parabolic velocity profile in regions 
where the boundary layer encompassed the entire film. The equation was 
discretized using the MacCormack predictor- corrector method [12]. For each 
time step, the predicted solution first used forward differencing in both 
time and space and then improved the prediction using forward differencing 
in time, but backward differencing in space. The corrected solution was 
taken as an arithmetic average of the past and predicted solutions. Since 
the one- dimensional solution did not require a significant amount of 
computer time, 500 grids were used in r- direction to obtain a smooth 
distribution of the film height. Convergence to the steady- state solution 
was obtained within 5000 iterations, where the maximum possible time- step 
that allowed stability was used. 

The complete two-dimensional solution of the flow and temperature 
fields were performed using a body- fitted curvilinear coordinate system. 
The grid cells were generated by algebraic interpolation between the 
boundaries of the computation domain, where the irregular free surface was 
taken as one of the boundaries. As shown in Fig. 1.1, the local 
coordinates were defined along lines connecting the centers of the adjacent 
grid cells. The discretized equations were written for each cell by using 
principles of conservation of mass, momentum and energy, where the primary 
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physical variables were retained without any nondimensionalization. The 
hybrid difference scheme presented by Patankar [13] was used to retain 
convection and/or diffusion terms. The flow field was solved using the 
SIMPLEST algorithm as demonstrated by Spalding [14] . The number of grid 
cells were determined from a series of tests that varied the number of 
grids in each coordinate direction. It was found that 50 x 30 cells in the 
z-y plane yielded grid independent solutions and was used for the present 
computation. 

The primary difficulty in two-dimensional computation was that the 
height of the free surface, which was unknown in general, needed to be 
specified in order to generate the grid structure and to solve the flow 
field. Therefore, an iterative procedure was required to ascertain the 
free surface height distribution. Two possible algorithms, outlined by 
Rahman et al. [9] and Rahman et al. [10] are ’pressure optimization method 
(POM)’ and ’porous wall method (PVM) ’ , respectively. Here, we would like 
to suggest improvements for both methods and give advantages of one over 
the other for certain kinds of problems. 

POM is suitable for one- domain flows where no jump is present, and 
two- domain flows where a jump isolates the supercritical and subcritical 
flow regimes. In this method, the free surface is approximated by an 
algebraic equation with a number of arbitrary constants. These constants 
are optimized by an exhaustive search technique which minimizes the error 
in the pressure distribution on the free surface defined as 
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Normalized RMS Error = 


1 n 
i E 

k = 1 


n 


Pg s. 


in 


+4 

Pk 

1 — ^ 
7 P V 


in 


where p^ is the difference between the free surface and ambient pressures 
at the kth radial location, and n is the number of cells in the radial 
direction. 


When a jump is present in the computation domain, the film height 
distribution is computed as follows: 


(1) The supercritical film height distribution for the entire disk is 
calculated by approximating the free surface by the equation 

6= S in ((r - A)/(r in - A)) 311 . 

The parameters ’A’ and ’an’ are optimized keeping the specified inlet 
film height constant. This equation allows 2 degrees of freedom. 
Higher degrees of freedom can be accomodated by assuming more complex 
forms of the equation. Note that supercritical flow does not require 
any specification of outlet condition, where the flow is assumed to be 
locally parabolic in nature. 

(2) The subcritical film height distribution for the entire disk is 
calculated by approximating the free surface by the equation 

• - S ont - »)/( r out - B » bn ' 

The parameters ’B’ and ’bn’ are optimized preserving the outlet film 
height fixed at determined by eqn. (1.11). A free falling 

condition was present at the outlet. 
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(3) The jump height distribution corresponding to the supercritical film 
height distribution is determined from the momentum balance at the 
jump [15] . 

if = 5 [ (1 * 9 -6 Fr?) 2 - l] 

where subscripts ’1’ and ’2’ indicate conditions upstream and 
downstream of the jump, respectively. 

(4) The location of the jump is determined by matching the jump height 
with corresponding subcritical film height distribution. 

(5) The supercritical film height distribution is recalculated for the 
partial disk from the entrance to the location of the jump. The 
method used is the same as that outlined in step (1) . 

(6) The slope of the jump is determined by calculating the flow as a 
single domain using the height determined in step (5) for the 
supercritical region and that from step (2) for the subcritical 
region. In this process, the jump surface is assumed to form a plane, 
the location for the initiation of jump is kept constant, and as 
before, the normalized RMS error is minimized. 

(7) The flow and energy equations are solved as a single domain problem 
using the film height distribution determined in step (6) to find the 
flow field and heat transfer coefficient. 

In steps (1), (2), (5) and (6), the film heights in the different 
regions are optimized independently to give a minimum normalized RMS error 
in pressure. The acceptable solution corresponds to a normalized RMS Error 


18 



of less than 0.1. This margin was established in [9] from the benchmark 
problem of falling film flow. Steps (5-7) are improvements over the 
original POM presented by Rahman et al. [9] , in which the jump occurred 
across one grid cell. So, the extent of the jump decreased with an 
increase in the number of computational cells. Moreover, the experiments 
of Thomas et al. [8] suggested a somewhat gradual jump rather than an 
instantaneous rise in the film height. Step (5) significantly improved the 
estimation of the film height in the supercritical region. Moreover, in 
step (2) , the free falling flow at the outlet was appropriately modeled and 
accounted for in the computation of subcritical film height distribution. 
In the original method [9] , the flow at the outlet was always assumed to be 
locally parabolic with a prescribed hydrostatic pressure gradient. In 
reality, the flow in a hydraulic jump is elliptic in nature, where outlet 
conditions significantly affect the location of the jump as well as the 
film height distribution. In order to accurately determine the 
heat- transfer coefficient in the vicinity of the jump, the one- domain 
calculation outlined in step (7) is essential since a natural continuity of 
the velocity and temperature fields needs to be preserved. It may also be 
noted that energy equation is solved only after the film height 
distribution is completely ascertained. 

In the PWM, the free surface is assumed to be a permeable wall through 
which fluid particles are allowed to enter or leave the control volume 
depending on the difference between the ambient pressure and the pressure 
in cells next to the free surface. The location of the surface is improved 
in successive iterations until the rate of penetration becomes negligible. 
Instead of assigning any particular equation for the free surface height 
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distribution, the new height is calculated from the old height and the 
local rate of penetration as shown below. 

. (Wk-WV f 1 - 12 *) 

(^corr)k = ( { old)k (1.12b) 

’in 

^newH ~ ^old^k + a ^corr^k (1.12c) 

where « is the relaxation factor. A value of a = 1 was found to be 
adequate for the problems considered here. The apparent advantage of this 
procedure is that there are more degrees of freedom for the free surface 
height distribution. The number of degrees of freedom here is the same as 
the number of grid cells in the radial direction. This technique is also 
more systematic and uniformly convergent than the exhaustive search 
technique associated with the POM. The deviation from the ideal 

zero- penetration condition can be estimated by either of the following 
measures . 


(a) Root- sum- square penetration 
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(b) Absolute sum of penetration 
^loss^kl 


n 

S I 
k = 1 


% 


m 


(c) Maximum error in flow rate 
I max 


(Q - »„) 


TJ 


in 


Since the method is uniformly convergent, all these error measures 
decrease with the number of iterations. In a practical problem, a solution 
may be acceptable when both root- sum- square penetration and maximum error 
in flow rate are less than 0.05 and the absolute sum of penetration is less 
than 0.1. 


The PVM outlined above is significantly different from the method 
originally proposed by Rahman et al. [10] . In the original method, the 
correction to the film height was obtained from the kinematic condition on 
the free surface instead of the local rate of penetration. This improved 
the surface in successive iterations to form a streamline, but apparently 
without a significant reduction in the local rate of penetration. 
Therefore, the convergence was extremely slow. Moreover, frequent 
smoothing was required to eliminate undesirable irregularities. These 
drawbacks were removed in the present modified algorithm. There is also a 
great improvement in the estimation of error. The criteria used by Rahman 
et al. [10] is RMS penetration, which is somewhat similar to the 
root- sum- square penetration except that an average was taken over all the 
cells in the radial direction. The disadvantage of using the average 
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penetration is that for a given free surface, it can be reduced by making 
the cell size smaller or increasing the number of cells in radial 
direction, which is somewhat misleading. In the present investigation, the 
total penetration over the entire free surface was considered, which was 
independent of the number of grid cells. Some limitations of the PVM were 
also discovered. It appears that the method does not work very well when a 
hydraulic jump is present, since a relatively large influx of fluid is 
encountered in the vicinity of the jump which leads to incorrect estimates 
of the subcritical film height downstream of the jump. However, it yields 
very good results when the flow is entirely supercritical, like flow in a 
zero- gravity environment . 

1.5 RESULTS AND DISCUSSION 

The film height distributions for thin film flows over a stationary 
disk are presented in Figures 1.2- 1.4 for flow rates of 7, 9, and 11 LPM, 
respectively. The inlet height corresponded to the first experimental data 
point, and the exit height was estimated by using equation (1.11). The 
two-dimensional solution used the modified POM, where the fluid velocity 
profile at the entrance was assumed to be parabolic with the maximum at the 
free surface. In the experiment by Thomas et al. [8], however, the fluid 
was discharged through a thin slot, so the maximum in the velocity profile 
at the discharge location is expected to be somewhere between the solid 
walls. After a short distance from the entrance, the flow is expected to 
evolve to a parabolic profile with the maximum velocity at the free 
surface, since the no- slip condition exists on the solid wall and the 
zero- shear condition is present on the free surface. Therefore, using the 
first measured data point instead of the physical entrance condition may be 
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Radial distance, mm 

Fig.1.2 Film height distribution for flow under normal gravity 

at a rate of 7 LPM 
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Radial distance, mm 

Fig.1.3 Film height distribution for flow under normal gravity 

at a rate of 9 LPM 



Radial distance, mm 

Flgl4 Film height distribution for flow under normal gravity 

at a rate of 11 LPM 




more appropriate for the numerical computations. 


At all flow rates, the present two-dimensional solution compares best 
with the experimental measurements. The present one- dimensional solution 
reasonably predicts the supercritical and subcritical film heights, but the 
jump appears to occur at a smaller radius. Moreover, the experimental data 
shows a somewhat gradual jump in contrast to the sudden jump predicted by 
the one- dimensional solution method. Figures 1.2- 1.4 also present a 
relative comparison of the present results with the one- dimensional 
solution of Thomas et al. [7] in which a critical condition (Fr=l) was 
specified at the outlet. The predicted subcritical film height 
distribution was much lower than that measured experimentally. The 
distribution of film height in the supercritical region is reasonably close 
to the experimental data and the two-dimensional solution. The solution of 
Thomas et al. [7] also predicted that the hydraulic jump occurred further 
downstream than the experiment. These comparisons show that the POM 
described above can be an effective mechanism for the prediction of the 
location of the jump as well as the supercritical and subcritical film 
height distributions. It may be improved further if the jump is assumed to 
be a curved surface instead of being a plane, as indicated by the data 
corresponding to 9 and 11 LPM. This, however, was not attempted here, 
since the shape of the film at the location of the jump is expected to be 
quite complicated due to churning motion of the fluid, and since no 
experimental data were available in regions where the jump is initiated. 
Figure 1.2 also shows a comparison with the analytical solution of Watson 
[1] for radially spreading thin film flow initiated by a liquid jet 
impinging perpendicularly to a solid wall. Since Watson’s results depend 
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on jet radius, an equivalent radius for the present system was calculated 


as a = v^2h c r c , where r £ is the radius and h £ is the thickness of the slot 
through which the liquid film is discharged. This equivalence was obtained 
from the relationship between jet radius and film thickness in an invisid 
flow. Vatson’s results show a larger film height than that measured 
experimentally or predicted numerically in the present investigation. The 
difference may be attributed to the approximate nature of the comparison 
due to the physical differences of the two flow systems. 


The details of the two-dimensional flow field in the presence of 
gravity is illustrated in Fig. 1.5 for the flow rate of 11 LPM. The 
velocity profile was parabolic at r = 0.101-0.110 m, which is in the 
supercritical region. Here the flow is well- structured with the maximum 
velocity at the free surface. This is the typical velocity profile in most 
of the supercritical region. However, as the jump is approached and film 
height increased rapidly, the inertia of the fluid stream decreased to 
preserve the continuity of fluid motion. Moreover, the stress- free 
condition on the free surface causes more forward- moving fluid particles to 
orient themselves to the free surface. Therefore, the flow separates from 
the solid wall, and a recirculating flow region is seen next to the wall 
just downstream of the jump. This is the characteristic seen in plots 
corresponding to r = 0.145-0.154 m and r = 0.174-0.183 m which are in the 
jump and suhcritical regions, respectively. The separation from the solid 
wall is due to frictional resistance in an adverse pressure gradient and 
was previously seen in the experimental work of Nakovyakov et al. [16]. 
They actually measured the wall shear stress for a circular hydraulic jump 
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Fig. 1.5 Velocity vectors for radially spreading flow under 
normal gravity at a rate of 11 I^M 



of a thin liquid film and found it to be negative for a short distance 
downstream from the location of the jump. 

The friction coefficient becomes zero at the point where the flow 
separates, and remains negative in the recirculating flow region. In the 
present work, it was seen that the flow remained separated in the entire 
subcritical region downstream of the jump. This is due to a relatively 
small extent of the subcritical region. A test run with a larger radius 
showed that the flow reattaches to the surface and a parabolic velocity 
profile is established further downstream. It should be also mentioned 
here that the above flow structure was obtained for a stationary free 
surface. The small variation of film height with time due to surface waves 
was neglected in the numerical computation. 

The variation of the film height for radially- spreading flow under 
zero gravity at a flow rate of 7 and 11 LPM are shown in Figures 1.6 and 
1.7, respectively. The dimension of the disk and the entrance flow 
condition were taken to be the same as the experiment of Thomas et al. [8] 
and numerical computation presented in Figs. 1.2 and 1.4. Since no 
experimental data is available for thin film flow in a gravity- free 
environment, the numerical prediction is the only means to understand the 
behavior of the flow. In the present work, the PVM was used to iterate the 
film height distribution, where the root- sum- square penetration and maximum 
excess flow rate were less than 5 percent, and absolute sum of penetration 
was less than 10 percent for the final solution. In a zero gravity 
environment, the film height increases monotonically with distance, and no 
sudden rise in the form of a hydraulic jump was encountered. Figures 1.6 
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Radial distance, mm 

Fig.1.6 Film height distribution for flow under zero gravity 

at a rate of 7 Ipm 








and 1.7 also show the height predicted by the one- dimensional computational 
method of Thomas et al. [7]. The one- dimensional solution gave a nearly 
linear rise in the film height, whereas the two-dimensional solution is 
curved. The difference between the two solutions is within 10 percent. It 
can be recalled here that the two-dimensional solution is obtained by 
successive reduction in the rate of penetration and the final film height 
distribution approached a streamline, which is an essential characteristic 
of the free surface. 

The distribution of velocity vectors for flow under zero gravity at a 
rate of 7 LPM is demonstrated in Fig. 1.8. The velocity profile across the 
thickness of the film was parabolic at all locations. The magnitude of the 
velocity, however, decreased as the fluid moved downstream. No jump or 
separation of flow was encountered in a zero-gravity environment, and the 
flow remained supercritical from entrance to exit. 

The distribution of Nusselt number for radially- spreading flow on a 
disk is shown in Fig. 1.9. The disk surface was maintained at a uniform 
temperature higher than the entering fluid temperature. The free surface 
was assumed to be adiabatic with no evaporation. The Nusselt number was 
calculated as Nu = htf/k, where h = q w /(T y - T^) . The Nusselt number 
decreased rapidly near the entrance due to the development of the thermal 
boundary layer from a uniform temperature condition at the entrance. The 
development of this boundary layer was extremely rapid. In the absence of 
gravity, the Nusselt number approached a uniform asymptotic value after a 
short distance from the entrance. In this case, the effect of flow rate 
was evident only in the thermally developing region and was very small in 
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Fig. 1.8 Velocity vectors for radially spreading flow under 
zero gravity at a rate of 7 LPM 
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Radial distance, mm 

Fig.1.9 Distribution of Nusselt number in radially spreading flow 






magnitude. The effect of flow rate became significant when gravity was 
present, since the hydraulic jump was encountered at different locations. 
In the supercritical region upstream from the jump, the Nusselt number 
approached a constant value following the development of thermal boundary 
layer. As soon as the jump was encountered, the Nusselt number first 
decreased, then increased to a peak, then decreased by a small amount, and 
increased monotonically further downstream. The irregular distribution in 
the vicinity of the jump is due to the churning motion of fluid and mixing 
of forward and recirculating flows. The monotonic increase in the 
subcritical region indicated the development of a thermal boundary layer in 
the recirculating flow region. The thickness of this boundary layer 
decreased with radius, and therefore, the Nusselt number increased. The 
maximum Nusselt number is seen to be at the exit end of the disk since the 
flow is forced to re- attach to the surface at that location due to the 
presence of the free falling liquid. Since the recirculating flow velocity 
was smaller than the main stream flow, the changes in the subcritical 
region were somewhat gradual in contrast to the rapid change in Nusselt 
number near the entrance. It can be recalled here that the film height, 
which was used as the lengthscale in the definition of the Nusselt number, 
was much larger in magnitude in the subcritical region than the 
supercritical region. To get a better understanding of the actual rate of 
heat transfer, the dimensional heat- transfer coefficient is presented in 
Fig. 1.10. 

The heat- transfer coefficient decreased monotonically downstream in a 
zero-gravity environment. In contrast to the Nusselt number, the 
heat- transfer coefficient was found to be a function of flow rate, and as 
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Fig.1.10 Distribution of local heat transfer coefficient 

in radially spreading flow 




expected, the rate of heat transfer increased with flow rate. In the 
presence of gravity, the heat- transfer coefficient decreased gradually 
downstream in the supercritical region, attained a sudden drop as soon as 
the jump was encountered, went through a somewhat irregular distribution 
and then increased gradually in the subcritical region following the jump. 
The lower heat- transfer coefficient in the vicinity of the jump was 
expected, since the main flow in the supercritical region and the 
recirculating flow in the subcritical region gather hot fluid there, 
leading to a smaller wall to fluid temperature gradient. The irregular 
variation is possibly the result of the churning motion of fluid particles 
in the vicinity of the jump. For the flow rates and dimensions considered 
here, the average heat- transfer coefficient over the entire plate was 
larger for the zero gravity flow than the normal gravity flow. Therefore, 
a better cooling of the disk may be expected when the system is operated in 
a zero or reduced gravity environment, as in a space shuttle. 

1.6 CONCLUSIONS 

The computational methods for the prediction of the film height and 
heat- transfer coefficient of free surface flow of a thin liquid film 
adjacent to a stationary disk were critically investigated, significantly 
improved, and compared with experimental measurements. A model for 
estimating the film height at the outer edge of the disk when gravitational 
body force is present was developed. The ’pressure optimization method’ 
was modified to appropriately handle flows involving a hydraulic jump. A 
new version of the ’porous wall method’ was presented, where iterations are 
more systematic, uniformly convergent, faster, and yield more accurate 
results. 
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In the absence of gravity, the film height increased monotonically , 
and the heat- transfer coefficient decreased gradually downstream. In this 
case, the Nusselt number approached a constant value after the initial 
development of thermal boundary layer. The velocity profile across the 
thickness of the film was parabolic with the maximum at the free surface. 

In the presence of gravity, a hydraulic jump was encountered that 
isolated the supercritical and subcritical flow regions. Flow separation 
was encountered at the location of the jump due to adverse pressure 
gradient and frictional resistance at the wall. The heat- transfer 
coefficient decreases downstream in the supercritical region and increases 
downstream in subcritical region with a sudden drop and fluctuations around 
the location of the jump. The predicted film height and location of the 
jump compared reasonably well with experimental measurements. 
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SECTION II 


FLOV OVER A ROTATING DISK 


2.1 SUMMARY 

The results of a numerical simulation of the flow field and associated 
heat transfer coefficient are presented for the free surface flow of a thin 
liquid film adjacent to a horizontal rotating disk. The computation has been 
performed for different flow rates and rotational velocities using a 
three-dimensional boundary- fitted coordinate system. Since the geometry of 
the free surface is unknown and dependent on flow rate, rate of rotation and 
other parameters, an iterative procedure had to be used to ascertain its 
location. The computed film height agreed well with existing experimental 
measurements. The flow is found to be dominated by inertia near the entrance 
and close to the free surface and dominated by centrifugal force at larger 
radii and adjacent to the disk. The rotation enhances the heat transfer 
coefficient by a significant amount. 
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2.2 INTRODUCTION 


The free surface flow of a thin liquid film adjacent to a rotating 
surface is an interesting fluid mechanics problem since a number of surface 
and body forces act on the system simultaneously to shape the flow structure. 
The most dominant of these forces are: the friction exerted by the disk; the 
centrifugal and Coriolis forces due to rotation of the disk; and the inertia 
from the incoming fluid stream. This kind of flow is also quite commonly 
encountered in engineering processes; e.g., evaporation or condensation on a 
turbine blade or spin coating of metals. The primary motivation of this 
study, however, is to understand the fluid flow and heat transfer in a 
proposed absorber unit for a space- based vapor absorption refrigeration 
system. This kind of system is expected to be very useful in a micro- gravity 
environment where the centrifugal body force can be an effective driving 
mechanism for the thinning of the film to promote the absorption of the 
refrigerant vapor into the absorbent. 

The fluid motion adjacent to a rotating surface in an infinite quiescent 
fluid medium is one of the fundamental problems in viscous fluid flow. 
Sparrow and Gregg [1] presented an exact solution of the Navier- Stokes 
equation for this case by estimating the thickness of the boundary layer from 
a basic mass and momentum balance and developing a similarity variable to 
transform the partial differential equations for the conservation of mass and 
momentum to a set of ordinary differential equations. These equations were 
integrated numerically across the thickness of the boundary layer to determine 
radial, axial and tangential velocities in dimensionless form. Murthy [2] 
analyzed the effects of turbulence in this flow by using the integral method 
and compared the results with available experimental data. 
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In the present study, however, we are not concerned with boundary layer 
flow in an infinite medium of fluid, but flow of a thin film adjacent to a 
plate. The boundary layer, however, develops inside this thin liquid film and 
may extend all the way across the film. An analytical solution for the flow 
of a thin film adjacent to a rotating disk was presented by Rauscher et al. 
[3] . An asymptotic expansion technique was used where the radial spread of 
fluid was perturbed to determine the effects of convection, Coriolis 
acceleration, radial diffusion, surface curvature and surface tension. These 
higher- order effects were discussed on a physical basis. Their solution was 
valid for laminar flow with small Rossby number. 

Espig and Hoyle [4] analytically studied the surface waves in a thin 
liquid layer on a rotating disk. Nusselt’s equation for film thickness on an 
inclined plate was modified to develop an expression for average film 
thickness, taking into account the effects of rotation. The thickness was 
found to be a function of volumetric flow rate, fluid viscosity and angular 
velocity of the disk. Needham and Merkin [5] studied the development of 
non-linear waves on the surface of a horizontally rotating thin liquid film. 
Using an asymptotic expansion technique, they studied the flow of a film much 
smaller in thickness than the distance from the axis of rotation. This 
distance, in turn, was smaller compared to the total radius of the disk in 
order to eliminate any end effect. It was found that the thickness of the 
film changes rapidly near the entrance due to spreading, while further away 
from the inlet, the film height is mainly determined by the centrifugal force 
and viscous stress. An expressions for film thickness was presented, which 
appears to be the same as that developed by Espig and Hoyle [4] . 
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The analysis of evaporation of a thin liquid film from a rotating surface 
was presented by Butuzov and Rifert [6]. A closed- form solution of the film 
thickness was derived by neglecting inertia and making a basic balance of 
centrifugal and frictional forces acting on the film. Bornside et al. [7] 
also studied evaporation from the free surface of a thin liquid film with 
reference to spin coating process. They developed a one- dimensional flow 
model accounting for variations of concentration, viscosity and diffusivity 
across the thickness of the spin coated film. The flow of the liquid was 
governed by a balance between centrifugal driving force and viscous resisting 
force. The equations of the flow and transport were solved using Galerkin 
finite element method. 

In all the above-mentioned studies, the thin film was primarily driven by 
centrifugal force. The inertia of the incoming fluid was either negligible or 
assumed to be so. In a space- based absorption refrigeration system, however, 
the film is expected to be introduced at the center of the disk with a 
significant amount of inertia to have a continuous flow of absorbent on the 
disk. This kind of system, where liquid is introduced at the center of the 
disk from a pressurized container and driven both by inertial and centrifugal 
forces, was studied recently by Thomas et al. [8,9]. 

Thomas et al. [8] presented the numerical solution of the one- 
dimensional radial flow of a thin liquid film adjacent to a stationary and 
rotating disk. The continuity and momentum equations were integrated across 
the thickness of the film by assuming a flat velocity profile to develop a 
single equation for film velocity. The resistance to the flow due to wall 
shear stress was expressed in terms of an empirical friction coefficient which 
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was estimated from the radial flow velocity. The thickness of the film was 
predicted by solving the equation for film velocity as a transient problem 
using the MacCormack predictor- corrector method. 

Thomas et al. [9] presented experimental measurements of radially 
spreading flow of a thin film adjacent to a horizontal disk. A non- obtrusive 
capacitance technique was used for the measurement of the film height 
distribution. A photographic study was done to understand the nature of the 
surface waves. Tests were performed both for stationary and rotating disks at 
flow rates ranging from 7 to 15 LPM. The rate of rotation varied between 0 
and 300 RPM. Their experiments predicted the presence of a hydraulic jump 
when the disk was stationary. The jump moved downstream with an increase of 
the flow rate and rotational velocity. Waves were found to be present on the 
free surface at all flow rates and rates of rotation. 

The present study is undertaken to develop a computational procedure for 
solving governing transport equations in their complete three-dimensional form 
for the flow of a thin liquid film adjacent to a rotating disk. This 
methodology elminated the need for assuming a friction coefficient as was done 
in the one- dimensional model developed by Thomas et al. [8]. Furthermore, an 
attempt is also made to improve the simple one- dimensional solution procedure. 
All computed results are compared with experimental measurements of Thomas et 
al. [10]. The details of the three-dimensional flow field and distribution of 
the heat- transfer coefficient are discussed. 

2.3 MATHEMATICAL MODEL 

The curvilinear boundary- fitted coordinate system used for the 
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three-dimensional numerical computation is shown in Fig. 2.1. The local 
coordinates are directed along lines connecting the centers of the adjacent 
grid cells. The x-axis is directed in the azimuthal direction, the y-axis 
perpendicular to the plate and the z-axis along the radial direction. The 
velocity components in these three directions are u, v and w, respectively. 
Due to the axisymmetric nature of the flow and vertical entrance and exit 
sections, u is directed parallel to the plate and v is directed perpendicular 
to the plate. However, the w- velocity changes its direction along the plate 
depending on the slope of the free surface. The height of the free surface 
from the solid wall is denoted by 6, which varies with radial location. 

For incompressible flow with constant fluid properties, the equations for 
the conservation of mass, momentum, and energy are given by 
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The terms due to viscous dissipation and pressure work in the energy 
equation are neglected as is typical for any low speed flow. Since the flow 
is symmetric about the axis of rotation, there was no variation of velocity or 
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temperature in the angular direction, i.e. d\/dx = 0 and dT/dx = 0. At the 
free surface, both tangential and normal stress components are zero. The 
vanishing tangential stress condition is represented by the zero velocity 
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Fig.2.1 Schematic of the flow of a thin liquid 
film adjacent to a rotating disk 
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gradient on the free surface. The normal stress condition on the free surface 
leads to equations balancing the pressure and other stresses, including 
surface tension. Scaling these equations, one can show that, for a thin film 
flow where Veber and Reynolds numbers are large, all other stress terms are 
negligible compared to pressure, leading to p = p Q on the free surface. The 
kinematic condition on the free surface relates the variation of film height 
to the velocity components. At the entrance plane, the velocity is assumed to 
be radial with a parabolic profile. At the exit plane, the flow is assumed to 
be fully developed with a hydrostatic pressure profile. The boundary 
conditions for the flow considered here are then given by 


at y = 0: 
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Here n is the direction normal to the free surface. The quantities T y , 

T T. , V- and a are assumed to be constant. Two different thermal 
sat ’ in’ in 

conditions are considered on the free surface: simple heating without 

evaporation; and evaporation. For the first case, the free surface is assumed 
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to be adiabatic in nature. In the second case, the free surface is isothermal 
at its equilibrium saturation temperature, T In the present 
investigation, it is also assumed that, in the case of evaporation, the fluid 
enters at its equilibrium saturation temperature, i.e., T^ n = T . . The 
Nusselt number is defined in terms of film thickness since that is the most 
significant length scale for this problem. 

2.4 COMPUTATIONAL PROCEDURE 

2.4.1 Three-Dimensional Method 

The governing transport equations along with the boundary conditions 
described in the previous section were solved numerically using a 
boundary- fitted curvilinear coordinate system. The irregular free surface was 
taken as one of the boundaries of the computational domain. The solution 
domain had a pie shape structure exteudiug from r in to r out iu the radial 
direction and over a small angle in the angular direction in order to prevent 
both distortion of body- fitted coordinates at large radii and clustering of 
grids at smaller radii. The grid cells were generated by an algebraic 
interpolation between the boundaries of the domain. In general, the cell 
faces were non- orthogonal to each other. As shown in Fig. 2.1, the local 
coordinates were defined along lines connecting the centers of the adjacent 
grid cells. The x-axis was taken in the azimuthal direction, the y-axis 
perpendicular to the plate and the z-axis in the radial flow direction. The 
velocity resolutes in these three directions are u, v, and w, respectively. 
In the finite- difference formulation, the co- variant components, i.e., 
components parallel to the cell faces, were used to represent the velocity and 
force vectors. 
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The finite- difference equations were derived by using the principle of 
conservation of mass, momentum and energy at each cell. The primary variables 
were preserved in the formulation instead of doing any non- dimens ionalizat ion. 
The quantities were stored in a staggered fashion where they made more 
physical sense for cell conservation. For each cell, the velocity components 
were stored at downstream boundaries, whereas all pressures and temperatures 
were stored at the cell center. The mass flux across a cell boundary was 
computed exactly from the scalar product of the velocity vector and the vector 
representing the area of the cell face. The convection contribution to a cell 
from its neighbors were calculated exactly by taking into account the 
curvature of a cell face and its non- orthogonal orientation. In the 
calculation of diffusion, however, the cell boundaries were approximated to be 
locally orthogonal. The hybrid difference scheme, demonstrated by Patankar 
[10] was used to preserve the relative contribution of convection and 
diffusion to a cell from its neighbor in terms of cell Peclet number. This is 
a common practice in the computation of convective flows. 

The distribution of cells in the computational domain was determined from 
a series of tests with different number of cells in the x, y, and z 
directions. Due to the axisymmetric nature of the flow, only 5 cells with an 
angular extent of 2.3°/cell were found to be adequate in that direction. This 
provided a total angle of 11.5° for the entire computation domain. Test runs 
with 3.45° /cell and 2.3° /cell yielded heat transfer coefficient within 0.001 
percent. Similarly, runs with 2 cells and 5 cells in the angular direction 
gave heat transfer coefficient within 0.0005 percent. Similar tests for the 
other two directions showed that 50 cells in the radial direction and 30 cells 
across the thickness of the film results in grid- independent solutions. 
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Therefore, the computation domain was divided into 5 x 30 x 50 cells by simple 
algebraic interpolation to generate a uniformly distributed grid structure. 

The finite difference equations were solved by using the SIMPLEST 

algorithm as presented by Spalding [11] , which is an iterative solution 

procedure where the computation was started by guessing a value for the 

pressure. Then the momentum equations were solved to determine the velocity 

components. The pressure was updated using the continuity equation. Even 

though the continuity equation does not contain any pressure, it can be easily 

transformed to a pressure correction equation as shown by Patankar [10] . The 

iterations were continued until the sum of the residuals for all computational 

- 

cells dropped below 10" for each equation. Since the flow was not coupled to 
the thermal transport for the problem considered here, the temperature field 
was solved once the film height distribution was completely ascertained. 

Since the geometry of the free surface was dependent on the flow 
parameters, but needed to be specified to generate the grid structure for the 
three-dimensional computation, an iterative procedure had to be used to 
ascertain the location of the free surface. The free surface boundary was 
assumed to be a permeable wall through which fluid particles were allowed to 
enter or leave the computational domain. Since the ambient pressure was 
prescribed, an outflow took place when the static pressure of the fluid 
adjacent to the free surface was higher than the ambient pressure and vice 
versa. The penetration of fluid through the free surface essentially violated 
the kinematic condition on the surface, which was arrived at by adjusting the 
surface height distribution in successive iterations. The scheme worked as 
follows: 
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1. A free surface height distribution was prescribed (based on 

one- dimensional solution of Thomas et al. [8] or any improved method. 

2. The flow field was solved completely for that height distribution. 

3. The amount of penetration of the fluid through the surface was calculated 

at all location along the flow. 

^loss^k = h-1 ' h 

4. The new free surface height was determined from the old height and the 
rate of penetration. 

^new^k = ^old^k + a Qin S k ^oldH 
where a is the relaxation factor. 

5. Iterations were continued until the rate of penetration became 

negligible. 

For a given surface height distribution, the deviation from the ideal 

zero penetration condition was estimated by the following measures. 

(a) Root square penetration 

k f x ^loss^k 
= «in 


(2.8a) 
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(b) Absolute sum of penetration 
n 



(c) Maximum error in flow rate 


Id’ 1in)U 



(2.8b) 


(2.8c) 


All these quantities were found to decrease almost monotonically with 
iterations. The final results presented here had both absolute sum of 
penetration and maximum error in flowrate of less than 0.06 and root square 
penetration less than 0.02. A relaxatim factor of a = 1 was found to be 
adequate in most computations. Since flow field was not affected by thermal 
transport, the energy equation was solved only for the final free surface 
height distribution with no penetration through the surface. The above 
computational algorithm is termed as porous wall method originally outlined by 
Rahman et al. [12]. Significant modifications of the method, however, have 
been made here to improve the iteration procedure as well as the error 
estimates. Moreover, this is the first time it has been applied to predict 
the three-dimensional flow involving rotation. 


2.4.2 Improved One- Dimensional Method 

The one- dimensional computational procedure, developed by Thomas et al. 
[8] was also improved here by incorporating a better estimate of frictional 
resistance exerted by the solid wall on the flow. In the one- dimensional 
method, the governing equations for mass and momentum were integrated across 
the thickness of the film assuming a uniform velocity profile. In the 
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original procedure of Thomas et al. [8], the flow was assumed to be strictly 
radial in nature with a superimposed solid- body rotation. The resistance to 
the flow due friction was taken to be the resistance in the radial direction. 
In the azimuthal direction, the velocity remained constant across the 
thickness of the film with no resistance form the solid wall. In reality, 
however, the velocity is expected to change due to the finite viscosity of the 
fluid, so there will be frictional resistance form the wall in the azimuthal 
direction. To account for the frictional resistance due to the angular 
velocity, we may define the total shear stress as 


r w 




+ r r 2 


(2.9) 


Here r f = c^ (1/2 p Ir), where c^ can be calculated by using the Blasius (see 
Kays and Crawford [13]) or parabolic solution given by Rahman et al. [14]. 
Also, can be estimated from the exact solution of laminar flow adjacent to 
a rotating disk in an infinite extent of fluid. As given by Schlichting [15], 


^ = 0.6 pr v l/2 # 3/2 (2.10) 

Using r w as the shear stress at the solid wall, the discretized equations of 
Thomas et al. [8] were reformulated and solved to give the film height 
distribution for any given flowrate and rate of rotation. 


2.5 RESULTS AND DISCUSSION 

The three-dimensional computational procedure discussed in the previous 
section was first tested against the limiting flow situations, where 
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theoretical or experimental results are already available. By setting the 
rotational speed to zero, the computed results reproduced exactly the film 
height distribution and velocity field corresponding to two-dimensional thin 
film flow adjacent to a stationary disk as presented by Rahman et al. [12]. 
Another limiting case was the flow motion induced by a rotating disk in an 
otherwise stagnant infinite fluid stream. Analytical solutions for this case 
was developed by Sparrow and Gregg [1]. In this system, the flow is limited 
to a thin boundary layer adjacent to the disk. The computed solution for this 
limiting case was plotted in Fig. 2.2, where it was compared with the 
analytical results of Sparrow and Gregg [1]. The non-dimensional distance, (, 
and velocity components G and F are defined as 


< = yJT 

(2.13) 

P _ W 

r o) 

(2.14) 

G = li- 
ra/ 

(2.15) 


It can been seen that both the radial and azimuthal velocity distributions 
across the thickness of the boundary layer are reasonably matched. One major 
difference between the conventional boundary layer flow and the flow of a thin 
film is that the thickness of the film is finite and bounded by the free 
surface which interacts with the surrounding medium. 

The height of a thin film adjacent to a rotating disk was experimentally 
measured by Thomas et al. [11]. A horizontal disk of diameter 406.4 mm was 
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used for the experiment where water at 20°C was introduced radially through a 
slot of height 0.267 mm at a radial location of 50.8 mm. The film height was 
measured by using a non-obstrusive capacitance probe from 76.2 mm to 195.6 mm. 
The measurements were not reported past 195.6 mm because part of the 
capacitance sensor was over the edge of the disk. Twenty sets of experimental 
data were taken covering flow rates of 7-15 LPM and rotational speeds of 
55-300 RPM. The flow was isothermal and no measurement of heat transfer was 
done in that experiment. In the present investigation, we used the flow and 
rotation conditions presented by Thomas et al. [11], so that a relative 
comparison of the film height distribution could be made. The specific cases 
chosen here are listed in Table 2.1. Prediction of film height and heat 
transfer coefficient could be done for all the experimental runs. However, to 
save the compuational efforts, some specific cases were picked up to cover the 
range of experiment in terms of flow rate and rate of rotation, as well as to 
understand the trend of variation of flow field and heat transfer coefficient 
with flow rate and rotational speed. The table also lists the root square 
penetration, absolute sum of penetration and maximum error in flow rate 
corresponding to the three-dimensional film height distribution for each 
specific flow rate and rate of rotation. For computational runs involving 
heat transfer, the wall temperature is assumed to be 30°C, whereas the 
entering fluid (water) is at 10°C. This gives a film temperature of 20°C for 
the evaluation of properties. 

The computed film height for a flow rate of 9 LPM and rotational speed of 
200 RPM (case 1 in Table 2.1) is shown in Fig. 2.3. The experimental data as 
well as the computed results using previous one- dimensional solution 
algorithms are also shown in the same figure. It can be noticed that the 
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Table 2.1 


Flow and Rotation Rates in Present Investigation 



Flow 

Rotational 


Absolute 

Maximum 


Rate 

Speed 

Root Square 

Sum of 

Error in 

Case & 

(LPM) 

(RPM) 

Penetration 

Penetration Flow Rate 

1 

9 

200 

0.0043 

0.0199 

0.0113 

2 

15 

200 

0.0072 

0.0235 

0.0152 

3 

7 

100 

0.0050 

0.0258 

0.0152 

4 

11 

100 

0.0110 

0.0406 

0.0218 

5 

15 

300 

0.0020 

0.0097 

0.0076 

6 

15 

55 

0.0127 

0.0593 

0.0540 

7 

7 

300 

0.0148 

0.0445 

0.0100 

8 

15 

100 

0.0092 

0.0332 

0.0224 
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Fig.23 Comparison of numerical solution with experimental data for 
thin liquid film at a flow rate of 9 LPM and angular velocity of 200 RPM 





three-dimensional solution predicts the trend of the experimental data, i.e., 
the film height first increases, attains a peak and then decreases further 
downstream. The increase of the film height near the entrance is due to 
strong frictional resistance the flow encounters from the solid wall. 
However, as the radius increases, the centrifugal force increases in magnitude 
and overpowers the frictional resistance. Then the film height decreases due 
to the increase in flow velocity as well as radial spreading. This latter 
effect is more significant at smaller radii. At larger radii, the curvature 
decreases and radially spreading flow approaches a plane flow. The 
one- dimensional solution of Thomas et al. [8] shows a monotonic decrease of 
film height with radius, and the height is consistently underpredicted except 
very close to the exit. Thomas et al. [8] assumed a solid- body rotation 
condition with no variation of velocity across the thickness of the film. 
Moreover, friction was computed based on the radial velocity using an 
empirical correlation. Therefore, the centrifugal force was more and 
frictional force was less accounted for in the numerical formulation. 
Therefore, the trend seen here is quite expected. The modified 
one- dimensional solution presented here improved the prediction of film height 
at smaller radii, but failed to do so at large radii. Here both radial and 
azimuthal components of frictional forces were accounted for, but still 
keeping no velocity variation across the thickness of the film. It may be 
noticed that the present three-dimensional solution, unlike the 
one- dimensional methods, does not require any estimate of friction 
coefficient. Figure 2.3 also shows a comparison with the analysis of Needham 
and Merkin [5] . Since the analysis presented an asymptotic solution valid for 
solid- body rotation, the analytical solution shows a monotic decrease in film 
height and compares well with the computation of Thomas et al. [8], where 
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similar assumptions were incorporated. In reality, however, viscous 
resistance at the wall slow down the film and result in a larger height 
distribution. Moreover, the asymptotic solution neglects the effects of 
inertia which was significant in the experiment and adequately taken care of 
in the three-dimensional solution. The computed results for a rotational rate 
of 200 RPM and flow rate of 15 LPM (case 2 in Table 2.1) are shown in Fig. 
2.4. Here, the present three-dimensional solution gives a very good agreement 
with the experimental data. The one- dimensional algorithm of Thomas et al. 
[8] consistently underpredicts the film height whereas the solution present 
one- dimensional algorithm improves it somewhat, particularly near the entrance 
and exit. The analytical solution of Needham and Merkin [5] shows a large 
film height near the entrance that monotonically decreases downstream. Near 
the entrance, the flow is dominated by inertia, and therefore, the analytical 
results of Needham and Merkin [5] is not valid in that region. 

The film height at a rotation rate of 100 RPM for two different flow 
rates (cases 3 and 4 in Table 2.1) are shown in Figs. 2.5 and 2.6. It can be 
seen that at this rate of rotation, the one- dimensional solution procedure 
using resultant frictional resistance gives a reasonable overall agreement, 
whereas the procedure of Thomas et al. [8] still underpredicts the film height 
at most locations. The analytical solution, as before, overpredicts the film 
height at smaller radii and underpredicts it at larger radii. The present 
three-dimensional solution gives the best agreement with the experimental data 
even though it does not coincide with it. A small overestimate of the film 
height near the entrance may be because of the development of the velocity 
profile from the given entrance boundary condition. At this boundary, the 
velocity was assumed to be parabolic in nature with a maximum on the free 
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Fig.24 Comparison of numerical solution with experimental data for 
thin liquid film at a flow rate of 15 LPM and angular velocity of 200 RPM 




Fig25 Comparison of numerical solution with experimental data for 
thin liquid film at a flow rate of 7 LPM and angular velocity of 100 RPM 







Rg2j6 Comparison of numerical solution with experimental data for 
thin liquid film at a flow rate of It LPM and angular velocity of 100 RPM 




surface. This is the best condition one can impose, particularly for laminar 
flow, as confirmed from the distribution of velocity across the thickness of 
the film in an intermediate location of the flow. 

The distribution of the film height for 300 RPM and flow rate of 15 LPM 
(case 5 in Table 2.1) is shown in Fig. 2.7, where it is compared with the 
one- dimensional predictions of Thomas et al. [8] experimental data of Thomas 
et al. [9] as well as the analysis of Needham and Merkin [5]. In this 
situation, the numerical solution underpredicts the experimental data in most 
regions. The one- dimensional algorithm of Thomas et al. [8] produces a film 
height distribution that is much smaller than the experimentally measured 
values. The present one- dimension algorithm improves it further, whereas, the 
present three-dimensional solution compares best with the experimental data. 
The trend of the data that the film height first increases, attains a peak and 
then decreases downstream is only captured by the present three-dimensional 
numerical solution. The one- dimensional solutions for this case show a 
monotonic decrease of film height with radius. It may be noticed here that 
analysis of Needham and Merkin [5] compares better with the computation of 
Thomas et al. [8] at a larger rate of rotation since the effects of inertia 
become smaller. The above comparisons (Figs. 2. 3-2. 7) show that the present 
algorithm can be an effective means for the prediction of the film height for 
a thin film flow adjacent to a rotating disk. The small deviations from the 
experimental values are within the uncertainties of experimental measurements 
and can also be attributed to unsteady surface waves as well as the surface 
tension effects that were not accounted for in the present steady- state 
numerical formulation, especially at the outer edge of the disk where the film 
height suddenly increased. 
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Fig27 Comparison of numerical solution with experimental data for 
thin liquid film at a flow rate of 15 LPM and angular velocity of 300 RPM 





The distribution of the velocity vectors at three locations across the 
film are shown in Figs. 2.8 and 2.9, respectively, for two limiting operating 
conditions: case 6 (55 RPM, 15 LPM) and case 7 (300 RPM, 7 LPM) . In Fig. 
2.8, it can be seen that the velocity was almost radial at y/6 = 0.88. Since 
the flow rate was high compared to the rate of rotation, the flow near the 
free surface was almost radial at this condition. At y/6 = 0.13, which was 
near the solid wall, it can be seen that the velocity vectors were directed 
more in the angular direction, particularly at larger radii. Looking at the 
magnitude of the velocity vectors, it can be seen that the centrifugal force 
increases as the radius increases. At y/6 = 0.48, it can be seen that the 
velocity vectors turn at an angle that increases with radial distance. Also, 
the overall magnitude of the vectors decreases with radius due to the 
resistance exerted by the solid wall. The flow is dominated by inertia at 
smaller radii and at locations away from the wall, whereas the contribution of 
rotation becomes important at locations near the solid wall and in regions 
away from the center. From a plot of the velocity components in the radial 
and azimuthal directions, it was found that the w- velocity profile is 
approximately parabolic in nature with a maximum at the free surface and zero 
at the wall. This corroborates corresponding observations for a stationary 
disk as presented by Rahman et al. [12]. The azimuthal component is maximum 
at the solid wall and gradually decreases as the free surface is approached. 
The normal component of velocity was found to be small compared to either of 
these components. 

In Fig. 2.9, it can be noticed that the flow is dominated by rotation in 
most regions except near the entrance when the rotational speed is 300 RPM and 
flow rate is 7 LPM. At this condition, the inertia is overpowered by the high 
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Fig.1-9 Velocity vectors for thin liquid film at a flow rate of 
7 LPM and angular velocity of 300 RPM 



rotational velocity after traversing some distance. Looking at the magnitude, 
it can also be noticed that the resultant velocity increases with radius. At 
intermediate conditions, the vector plots showed behavior within these two 
limits. The effect of the flow rate was found to increase the radial 
component of velocity, whereas rotation increased both radial and azimuthal 
components. The friction always counteracts these velocities to arrive at an 
equilibrium. The fluid particles exit the disk at an angle with the radial 
direction that increases with the rate of rotation. 

A plot of Nusselt number for a flow rate of 15 LPM at two different rates 
of rotation are shown in Fig. 2.10 (cases 5 and 8 in Table 2.1). The plots 
correspond to two different thermal boundary conditions on the free surface, 
namely, simple heating without evaporation, when the free surface is 
approximately adiabatic in nature, and evaporative free surface when the 
surface is isothermal at the equilibrium saturation temperature. The thermal 
condition at the rotating disk surface is assumed to be a uniform temperature 
higher than the entering fluid (or saturation) temperature. It can be seen 
that the Nusselt number decreases rapidly near the entrance at all flow and 
heating conditions. This is due to the development of the thermal boundary 
layer as the fluid moves downstream. For both the cases of heating and 
evaporation, it can be seen that the Nusselt number for the case of 300 RPM 
decreases monotonically all the way to the exit, whereas in the case of 100 
RPM, the Nusselt number decreases first, attains a minimum and then increases 
further downstream to approach an approximately constant value near the exit. 
The Nusselt number for heating is consistently found to be higher than the 
corresponding case of evaporation at all radial locations and all rates of 
rotation. This is also observed in connection with heat transfer to a falling 
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Fig.210 Nusselt number distribution for thin liquid film at 

a flow rate of 15 LPM 





liquid film presented by Seban and Faghri [16] . The Nusselt number for 300 
RPM is found to be higher than that for 100 RPM in most regions of the disk 
except close to the exit. The slightly higher Nusselt number near the exit is 
due to much larger film height than the corresponding case of 300 RPM. From 
this figure it can also be seen that an overall enhancement of Nusselt number 
is obtained by increasing the rate of rotation. 

The Nusselt number here is defined in terms of local film thickness, 
which also varies along the radius. The film thickness also decreases as the 
rotational velocity increases. To see the distribution of the heat- transfer 
rate more clearly, the values of the heat- transfer coefficient for the same 
flow rate and heating conditions are shown in Fig. 2.11. It can be seen that 
heat- transfer coefficient for 300 RPM is consistently higher than that for 
100 RPM at all radial locations. This is quite expected since the rate of 
rotation increases the centrifugal force which accelerates and thins the film, 
both of which contribute to the enhancement of heat transfer. It can also be 
seen that for the case of 300 RPM and simple heating, the heat transfer 
coefficient decreases near the entrance, attains a minimum and then increases 
further downstream. For the other cases shown here, the heat transfer 
coefficient approaches approximately a constant value after gradual decrease 
near the entrance. As the radius increases, the inertial force decreases and 
centrifugal force increases, both of which are counteracted by frictional 
resistance at the solid wall. In the absence of rotation or any external body 
force, the heat- transfer coefficient gradually decreases downstream as seen in 
the studies by Rahman et al. [12] for radial flow at zero gravity. The 
rotation balances out or even overcomes the frictional resistance at larger 
rate of rotation, as seen in the vector plot for 300 RPM (Fig. 2.9). 
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Therefore, the convective heat transfer, which is intimately related to the 
fluid velocity, increases downstream once the centrifugal force becomes the 
dominant driving mechanism. If a disk of larger radius was considered, we 
could possibly see a minimum followed by a gradual rise in the heat- transfer 
coefficient for all cases. It can also be noticed that the effects of 
rotation in the enhancement of the heat- transfer coefficient is more for the 
case of simple heating than the case of evaporation. 

Figure 2.12 shows the variation of Nusselt number with radial distance 
for a given rate of rotation at different flow rates. For the case of 15 LPM, 
the Nusselt number monotonically decreases from the entrance radius, but for 
the case of 7 LPM, it decreases gradually at smaller radii and approaches a 
constant value further downstream. The behavior is the same for both the 
cases of heating and evaporation. For the disk considered here, the Nusselt 
number for 15 LPM is found to be higher than that for 7 LPM at all radial 
locations. This is expected, particularly at smaller radial locations, since 
fluid velocity will be higher at larger flow rates. At large radii, however, 
when effects of rotation become more important, a higher heat- transfer 
coefficient may be attained by using a smaller flow rate due to more thinning 
of the film. So, by controlling the disk radius, flow rate and rate of 
rotation, one may attain any desired requirement of heat transfer in this kind 
of flow system. 

2.6 CONCLUSIONS 

Numerically computed distributions of the film height, velocity vectors 
and heat- transfer coefficient are presented for the free surface flow of a 
thin liquid film adjacent to a horizontal rotating disk. The flow is found to 
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be dominated by inertia at smaller radii and close to the free surface and 
dominated by rotation near the solid wall and at larger radii. The radial 
component of velocity is found to have an approximately parabolic profile 
which is maximum at the free surface and zero at the solid wall. The angular 
component of velocity is maximum at the wall and gradually diminishes across 
the thickness of the film. The local fluid velocity is increased with an 
increase in both the volumetric flow rate and angular velocity of rotation. 
The fluid particles exit the disk at an angle with radius that increases with 
rate of rotation. At the flow rates and rates of rotation considered here, 
the film height first increases, attains a peak and then decreases further 
downstream. The increment of the film height is attributed to frictional 
resistance whereas the reduction of the film height at larger radii is due to 
the spreading of the film as the flow area increases as well as the increase 
of the centrifugal force with radius. The predicted height using the present 
three-dimensional computational method agreed well with experimental 
measurements. The rate of heat transfer is enhanced by a significant amount 
by increasing the rate of rotation for both the cases of simple heating with 
no evaporation and evaporation at the free surface. 
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